home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / voronoi.pro < prev    next >
Text File  |  1997-07-08  |  8KB  |  247 lines
  1. ; $Id: voronoi.pro,v 1.7 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1992-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. function isright, x0, y0, x1, y1, x2, y2 
  7. ; return 1 if Pnt0 is to right of Pnt1-> Pnt2
  8. ; return 0 if it is on the line.
  9. ; return -1 if Pnt0 is to the left of Pnt1 -> Pnt2
  10.  
  11. z = (x0-x1) * (y2-y1) - (y0-y1) * (x2-x1)
  12. if z gt 0.0 then return, 1
  13. if z lt 0.0 then return, -1
  14. return, 0
  15. end
  16.  
  17.  
  18. PRO voronoi_get_intersect, rect, x0, y0, dx, dy, xi, yi, nedge
  19. ;  Return the closest intersection of the line thru (x0, y0), with
  20. ;    slope (dx / dy), and the rectangle rect.  The intersection must
  21. ;    be in the "positive" direction.
  22. ;    Set (xi, yi) to the intersection point, and 
  23. ;    nedge to the edge index.  Edge 0 is the bottom, going CCW.
  24. ;
  25. nedge = -1
  26. tmin = 0.
  27.  
  28. if dy ne 0.0 then begin
  29.     t = (rect[1] - y0) / dy    ;Bottom = edge 0?
  30.     if t ge 0.0 then begin
  31.     tmin = t
  32.     nedge = 0
  33.     yi = rect[1]        ;Intersection with bottom
  34.     xi = x0 + t * dx
  35.     endif
  36.     t = (rect[3] - y0) / dy    ;Top edge = 2
  37.     if t ge 0.0 then $
  38.       if (nedge lt 0) or (t lt tmin) then begin
  39.     tmin = t
  40.     nedge = 2
  41.     yi = rect[3]
  42.     xi = x0 + t * dx
  43.     endif
  44. endif                ;dy ne 0
  45.  
  46. if dx ne 0.0 then begin        ;Check sides?
  47.     t = (rect[0] - x0) / dx    ;Left edge = 3
  48.     if t ge 0.0 then $
  49.        if (nedge lt 0) or (t lt tmin) then begin
  50.     tmin = t
  51.     nedge = 3
  52.     xi = rect[0]
  53.     yi = y0 + t * dy
  54.     endif
  55.     t = (rect[2] - x0) / dx    ;Right edge = 1
  56.     if t ge 0.0 then $
  57.        if (nedge lt 0) or (t lt tmin) then begin
  58.     tmin = t
  59.     nedge = 1
  60.     xi = rect[2]
  61.     yi = y0 + t * dy
  62.     endif
  63. endif                ;Dx ne 0
  64. end
  65.  
  66.  
  67.  
  68.  
  69. pro VORONOI_SHOW, n        ;Illustrate using and drawing Voronoi polygons
  70. ; This procedure generates N random points (default = 12).
  71. ;    and then draws the voronoi polygons, with the points and
  72. ;    delaunay triangulation overlaid.
  73.  
  74. if n_elements(n) le 0 then n = 12    ;Make the random points
  75. seed = 1211567L
  76. x = randomu(seed, n)
  77. y = randomu(seed, n)
  78.  
  79. triangulate, x, y, tr, CONN=c        ;Triangulate them
  80. plot,x,y,/psym, xrange=[-.5,1.5], yrange=[-.5,1.5]
  81. tek_color                ;Discrete color tables.
  82. range = 0                ;Init bounding rectangle
  83.  
  84. for i=0, n-1 do begin        ;Each VORONOI region for each point
  85.     voronoi, x, y, i, c, xp, yp, range      ;Get the ith polygon
  86.     xp = xp > (-10) < 10        ;Clip to reasonable space
  87.     yp = yp > (-10) < 10
  88.     polyfill, xp, yp, color = (i mod 13) + 2    ;Show it
  89.     endfor
  90.  
  91. oplot,x,y,/psym, color=1        ;Show points & triangulation
  92. for i=0, n_elements(tr)/3-1 do begin    ;The triangles
  93.     t = [tr(*,i),tr(0,i)]        ;Subscripts of triangle & back to 0
  94.     plots,x(t), y(t)
  95.     endfor
  96. end
  97.  
  98.  
  99.  
  100. PRO voronoi, x, y, i0, c, xp, yp, rect
  101. ; Copyright (c) 1992, Research Systems, Inc. All rights reserved.
  102. ;    Unauthorized reproduction prohibited.
  103. ;+
  104. ; NAME:
  105. ;    VORONOI
  106. ;
  107. ; PURPOSE:
  108. ;    This procedure computes the Voronoi polygon of a point within
  109. ;    an irregular grid of points, given the Delaunay triangulation.
  110. ;    The Voronoi polygon of a point contains the region closer to
  111. ;    that point than to any other point.
  112. ;
  113. ; CATEGORY:
  114. ;    Gridding.
  115. ;
  116. ; CALLING SEQUENCE:
  117. ;    VORONOI, X, Y, I0, C, Xp, Yp, Rect
  118. ;
  119. ; INPUTS:
  120. ;    X:    An array containing the X locations of the points.
  121. ;    Y:    An array containing the Y locations of the points.
  122. ;    I0:    Index of the point of which to obtain the Voronoi polygon.
  123. ;    C:    A connectivity list from the Delaunay triangulation.
  124. ;        This list is produced with the CONNECTIVITY keyword
  125. ;        of the TRIANGULATE procedure.
  126. ;    Rect    the bounding rectangle:  [Xmin, Ymin, Xmax, Ymax].
  127. ;        Because the Voronoi polygon (VP) for points on the convex hull
  128. ;        extends to infinity, a clipping rectangle must be supplied to
  129. ;        close the polygon.  This rectangle has no effect on the VP of
  130. ;        interior points.  If this rectangle does not enclose all the
  131. ;        Voronoi vertices, the results will be incorrect.  If this
  132. ;        parameter, which must be a named variable, is undefined or
  133. ;        set to a scalar value, it will be calculated.
  134. ;
  135. ; OUTPUTS:
  136. ;    Xp, Yp:    The vertices of voroni polygon, VP.
  137. ;
  138. ; RESTRICTIONS:
  139. ;    The polygons only cover the convex hull of the set of points.
  140. ;
  141. ; PROCEDURE:
  142. ;    For interior points, the polygon is constructed by connecting
  143. ;    the midpoints of the lines connecting the point with its Delaunay
  144. ;    neighbors. Polygons are traversed in a counterclockwise direction.
  145. ;
  146. ;    For exterior points, the set described by the midpoints of the
  147. ;    connecting lines, plus the circumcenters of the two triangles
  148. ;    that connect the point to the two adjacent exterior points.
  149. ;
  150. ; EXAMPLE:
  151. ;    See the example procedure, VORONOI_SHOW, contained in this file.
  152. ;    To illustrate Voronoi polygons, after compiling this file (voronoi):
  153. ;        VORONOI_SHOW, Npoints  (try anywhere from 3 to 1000, default=12)
  154. ;
  155. ;    To draw the voroni polygons of each point of an irregular 
  156. ;    grid:
  157. ;      x = randomu(seed, n)             ;Random grid of N points
  158. ;      y = randomu(seed, n)
  159. ;      triangulate, x, y, tr, CONN=c              ;Triangulate it
  160. ;      rect = 0
  161. ;      for i=0, n-1 do begin
  162. ;        voronoi, x, y, i, c, xp, yp, rect      ;Get the ith polygon
  163. ;        polyfill, xp, yp, color = (i mod 10) + 2  ;Draw it
  164. ;        endfor
  165. ;
  166. ; MODIFICATION HISTORY:
  167. ;    DMS, RSI.    Dec, 1992. Original version.
  168. ;    DMS, RSI    Feb, 1995. Added bounding rectangle which simplified
  169. ;                   logic and better illustrated VPs for points
  170. ;                   on the convex hull.
  171. ;-
  172. COMMON VORONOI_COMMON, first    ;Only print warning once.
  173.  
  174. if n_params() lt 7 and n_elements(first) le 0 and !quiet eq 0 then begin
  175.   Message, /INFO, 'New revision. For more efficient use, supply the Rect parameter'
  176.   first = 1
  177.   endif
  178.  
  179. p = c(c(i0):c(i0+1)-1)      ;Verts of polygon
  180. np = n_elements(p)
  181. xp = fltarr(np, /NOZERO)
  182. yp = fltarr(np, /NOZERO)
  183. ext = i0 eq p(0)        ;True if exterior point
  184.  
  185. ; Each vertex is simply the circumcenter of a Delaunay triangle
  186. ;         containing the point in question.
  187. for i=0, np-1 do begin        ;Traverse adjacency list for point i0.
  188.     m = p(i)
  189.     j = p((i + 1) mod np)    ;Successor
  190.     cir_3pnt, x([m,j,i0]), y([m,j,i0]), r, x0, y0
  191.     xp[i] = x0
  192.     yp[i] = y0
  193.     endfor
  194.  
  195. if ext eq 0 then return    ;If interior point, we're all done....
  196.  
  197. ;    *** Point is on the Convex Hull ****
  198. if n_elements(rect) ne 4 then begin  ;Initialize bounding rect?
  199.     i1 = i0            ;Follow boundary CCW
  200.     xmin = min(x, max=xmax)
  201.     ymin = min(y, max=ymax)
  202.     xr = (xmax-xmin)*.05    ;Fudge factor of 5%
  203.     yr = (ymax-ymin)*.05
  204.     rect = [xmin-xr, ymin-yr, xmax+xr, ymax+yr]  ;Initial bound box
  205.     repeat begin        ;Get circumctr of each boundary triangle
  206.     k = c(i1)
  207.     m = c(k+1)
  208.     j = c(k+2)
  209.     cir_3pnt, x([m,j,i1]), y([m,j,i1]), r, x0, y0
  210.     rect[0] = rect[0] < x0
  211.     rect[2] = rect[2] > x0
  212.     rect[1] = rect[1] < y0
  213.     rect[3] = rect[3] > y0    ;Keep extremes
  214.     i1 = m            ;Next boundary point
  215.     endrep until i1 eq i0
  216.     endif
  217.     
  218. ;    Now get intersections of perpendicular  bisectors of the two edges
  219. ;    on the convex hull, with vertex point i0, with the bounding rectangle.
  220. ;
  221. j = p(1)
  222. voronoi_get_intersect, rect, xp[1], yp[1], y[j]-y[i0], x[i0]-x[j], x0, y0,edge0
  223. xp[0] = x0
  224. yp[0] = y0
  225. j = p(np - 1)
  226. voronoi_get_intersect, rect, xp[np-2], yp[np-2], y[i0]-y[j], x[j]-x[i0], $
  227.     x0, y0, edge1
  228. xp[np-1] = x0
  229. yp[np-1] = y0
  230.  
  231. if (edge0 < edge1) lt 0 then begin   ;Either out of bounds?
  232.     MESSAGE, /INFO, 'Bounding rectangle does not enclose Voronoi polygon'
  233.     xp[0] = xp[1]        ;Fudge polygon, its wrong anyway...
  234.     yp[0] = yp[1]
  235.     xp[np-1] = xp[np-2]
  236.     yp[np-1] = yp[np-2]
  237.     return
  238. endif
  239.  
  240. while edge1 ne edge0 do begin        ;Add corner(s) of bound rect if necess.
  241.     edge1 = (edge1 + 1) mod 4        ;Go CCW
  242.     xp = [xp, rect[([0,2,2,0])[edge1]]]
  243.     yp = [yp, rect[([1,1,3,3])[edge1]]]
  244.     endwhile
  245. return
  246. end
  247.